Ejemplo de análisis espacial: densidad de la red vial#
Introducción#
La densidad de la red vial para un polígono se define como:
\[longitudRedVial(km) / áreaPolígono(km^{2})\]
Por ejemplo, si un cantón tiene 500 km de longitud de red vial y un área de 1000 km2, la densidad de su red vial es 0.5.
Este cuaderno de notas calcula la densidad de la red vial para cada uno de los cantones de Costa Rica y la presenta en un mapa de coropletas.
Carga de bibliotecas#
# Carga de geopandas con el alias gdp
import geopandas as gpd
# Otras bibliotecas
# Carga de pandas con el alias pd
import pandas as pd
# Carga del módulo pyplot de matplotlib con el alias plt
import matplotlib.pyplot as plt
# Carga de la clase WebFeatureService del módulo wfs de owslib
# Permite interactuar con servicios web geoespaciales tipo WFS
from owslib.wfs import WebFeatureService
# Carga de la clase BytesIO del módulo estándar io
# Permite crear un objeto en memoria que actúa como un archivo binario
from io import BytesIO
# Carga de la biblioteca Folium, para mapas interactivos
import folium
Carga de datos#
Cantones#
# Conexión al servicio WFS del IGN 5k CO
wfs_url = 'https://geos.snitcr.go.cr/be/IGN_5_CO/wfs'
wfs_version = '1.1.0'
wfs = WebFeatureService(url=wfs_url, version=wfs_version)
# Obtener la capa de cantones
capa = 'IGN_5_CO:limitecantonal_5k'
respuesta = wfs.getfeature(typename=capa, outputFormat='application/json')
# Leer la respuesta en un GeoDataFrame
cantones_gdf = gpd.read_file(BytesIO(respuesta.read()))
# Reducción de columnas
cantones_gdf =cantones_gdf[['CÓDIGO_CANTÓN', 'CANTÓN', 'geometry']]
# Definir un índice
# cantones_gdf.set_index('CÓDIGO_CANTÓN', inplace=True)
# Mostrar las primeras filas del GeoDataFrame
print(cantones_gdf.head())
# Mapa
cantones_gdf.plot()
CÓDIGO_CANTÓN CANTÓN \
0 101 San José
1 102 Escazú
2 103 Desamparados
3 104 Puriscal
4 105 Tarrazú
geometry
0 MULTIPOLYGON (((481030.328 1102633.167, 481032...
1 MULTIPOLYGON (((479844.919 1102669.584, 479851...
2 MULTIPOLYGON (((494091.231 1094936.82, 494091....
3 MULTIPOLYGON (((455934.039 1095770.324, 455949...
4 MULTIPOLYGON (((501998.931 1074558.189, 502000...
<Axes: >
Red vial#
# Conexión al servicio WFS del IGN 200 k
wfs_url = 'https://geos.snitcr.go.cr/be/IGN_200/wfs'
wfs_version = '1.1.0'
wfs = WebFeatureService(url=wfs_url, version=wfs_version)
# Obtener la capa de red vial
capa = 'IGN_200:redvial_200k'
respuesta = wfs.getfeature(typename=capa, outputFormat='application/json')
# Leer la respuesta en un GeoDataFrame
red_vial_gdf = gpd.read_file(BytesIO(respuesta.read()))
# Añadir una columna con un identificador único
red_vial_gdf['red_vial_id'] = range(1, len(red_vial_gdf) + 1)
# Reducción de columnas
red_vial_gdf =red_vial_gdf[['red_vial_id', 'geometry']]
# Definir un índice
# red_vial_gdf.set_index('red_vial_id', inplace=True)
# Mostrar las primeras filas del GeoDataFrame
print(red_vial_gdf.head())
# Mapa
red_vial_gdf.plot()
/home/mfvargas/miniconda3/envs/gf0657-2024-ii/lib/python3.12/site-packages/pyogrio/raw.py:196: RuntimeWarning: Several features with id = 0 have been found. Altering it to be unique. This warning will not be emitted anymore for this layer
return ogr_read(
red_vial_id geometry
0 1 LINESTRING (543973.13 1157346.96, 543972.52 11...
1 2 LINESTRING (542005.8 1160843.73, 542304.68 116...
2 3 LINESTRING (542872.85 1160705.42, 542946.13 11...
3 4 LINESTRING (527814.53 1167099.79, 527692.44 11...
4 5 LINESTRING (538428.65 1157738.05, 538251.35 11...
<Axes: >
Análisis#
# Generar un segmento de red vial
# por cada intersección con un polígono de cantón
interseccion_redvial_cantones_gdf = gpd.overlay(
red_vial_gdf,
cantones_gdf,
how='intersection'
)
interseccion_redvial_cantones_gdf.head()
/tmp/ipykernel_228903/3535146563.py:3: UserWarning: CRS mismatch between the CRS of left geometries and the CRS of right geometries.
Use `to_crs()` to reproject one of the input geometries to match the CRS of the other.
Left CRS: EPSG:5367
Right CRS: EPSG:8908
interseccion_redvial_cantones_gdf = gpd.overlay(
| red_vial_id | CÓDIGO_CANTÓN | CANTÓN | geometry | |
|---|---|---|---|---|
| 0 | 1 | 702 | Pococí | LINESTRING (543973.13 1157346.96, 543972.52 11... |
| 1 | 2 | 702 | Pococí | LINESTRING (542005.8 1160843.73, 542304.68 116... |
| 2 | 3 | 702 | Pococí | LINESTRING (542872.85 1160705.42, 542946.13 11... |
| 3 | 4 | 702 | Pococí | LINESTRING (527814.53 1167099.79, 527692.44 11... |
| 4 | 5 | 702 | Pococí | LINESTRING (538428.65 1157738.05, 538251.35 11... |
# Calcular la longitud de cada segmento de red vial intersecado
interseccion_redvial_cantones_gdf["LONGITUD_RED_VIAL"] = interseccion_redvial_cantones_gdf["geometry"].length/1000
interseccion_redvial_cantones_gdf.head()
| red_vial_id | CÓDIGO_CANTÓN | CANTÓN | geometry | LONGITUD_RED_VIAL | |
|---|---|---|---|---|---|
| 0 | 1 | 702 | Pococí | LINESTRING (543973.13 1157346.96, 543972.52 11... | 8.744950 |
| 1 | 2 | 702 | Pococí | LINESTRING (542005.8 1160843.73, 542304.68 116... | 8.188802 |
| 2 | 3 | 702 | Pococí | LINESTRING (542872.85 1160705.42, 542946.13 11... | 2.135535 |
| 3 | 4 | 702 | Pococí | LINESTRING (527814.53 1167099.79, 527692.44 11... | 2.302054 |
| 4 | 5 | 702 | Pococí | LINESTRING (538428.65 1157738.05, 538251.35 11... | 8.779122 |
# Agrupación y suma de los los segmentos de red vial por cantón
suma_redvial_cantones_df = interseccion_redvial_cantones_gdf.groupby(["CÓDIGO_CANTÓN"])["LONGITUD_RED_VIAL"].sum()
# Convertir la serie en dataframe
suma_redvial_cantones_df = suma_redvial_cantones_df.reset_index()
# Desplegar longitudes de red vial por cantón
suma_redvial_cantones_df.sort_values(by='LONGITUD_RED_VIAL', ascending=False)
| CÓDIGO_CANTÓN | LONGITUD_RED_VIAL | |
|---|---|---|
| 29 | 210 | 2591.051204 |
| 18 | 119 | 1429.548737 |
| 65 | 601 | 1319.628280 |
| 67 | 603 | 1297.054690 |
| 53 | 410 | 1103.343280 |
| ... | ... | ... |
| 14 | 115 | 23.190868 |
| 50 | 407 | 20.224896 |
| 12 | 113 | 16.403583 |
| 52 | 409 | 15.086709 |
| 51 | 408 | 11.006955 |
84 rows × 2 columns
# Agregar a la capa de cantones la longitud de la red vial para cada polígono
cantones_redvial_gdf = cantones_gdf.merge(suma_redvial_cantones_df, on="CÓDIGO_CANTÓN")
# Calcular el área del cantón
cantones_redvial_gdf["ÁREA_CANTÓN"] = cantones_redvial_gdf["geometry"].area/1000000
# Calcular la densidad de la red vial
cantones_redvial_gdf["DENSIDAD_RED_VIAL"] = cantones_redvial_gdf["LONGITUD_RED_VIAL"] / cantones_redvial_gdf["ÁREA_CANTÓN"]
# Desplegar la densidad de red vial en cada cantón
cantones_redvial_gdf[['CÓDIGO_CANTÓN', 'CANTÓN', 'DENSIDAD_RED_VIAL', 'ÁREA_CANTÓN', 'LONGITUD_RED_VIAL']].sort_values(by="DENSIDAD_RED_VIAL", ascending=False)
| CÓDIGO_CANTÓN | CANTÓN | DENSIDAD_RED_VIAL | ÁREA_CANTÓN | LONGITUD_RED_VIAL | |
|---|---|---|---|---|---|
| 0 | 101 | San José | 2.259410 | 44.624284 | 100.824537 |
| 7 | 108 | Goicoechea | 2.148213 | 31.696493 | 68.090826 |
| 17 | 118 | Curridabat | 2.065924 | 16.062606 | 33.184129 |
| 12 | 113 | Tibás | 1.983644 | 8.269421 | 16.403583 |
| 19 | 120 | León Cortés Castro | 1.954404 | 121.894583 | 238.231301 |
| ... | ... | ... | ... | ... | ... |
| 69 | 605 | Osa | 0.351286 | 1932.697429 | 678.930001 |
| 78 | 701 | Limón | 0.291557 | 1769.378451 | 515.874969 |
| 77 | 613 | Puerto Jiménez | 0.287517 | 720.432185 | 207.136511 |
| 44 | 401 | Heredia | 0.287120 | 283.109137 | 81.286423 |
| 81 | 704 | Talamanca | 0.117583 | 2792.225947 | 328.318153 |
84 rows × 5 columns
# Mapa estático
cantones_redvial_gdf.plot(
column="DENSIDAD_RED_VIAL",
legend=True,
cmap='OrRd',
scheme='quantiles',
figsize=(20, 20)
)
<Axes: >
# Antes de agregar cantones_redvial_gdf,
# se simplifica (i.e. se reduce el número de vértices)
# para generar un HTML más pequeño (< 100 MB)
# Hacer una copia de cantones_redvial_gdf para simplificarla
cantones_redvial_simplificado_gdf = cantones_redvial_gdf
# Definir el factor de simplificación
# Cuanto mayor sea el valor, mayor será la simplificación
factor_simplificacion = 0.01
# Simplificar geometría
cantones_redvial_simplificado_gdf['geometry'] = cantones_redvial_simplificado_gdf['geometry'].simplify(
factor_simplificacion,
preserve_topology=True
)
# Crear el mapa interactivo con la densidad de la red vial
m = cantones_redvial_simplificado_gdf.explore(
column='DENSIDAD_RED_VIAL',
name='Densidad de la red vial en cantones',
cmap='OrRd',
tooltip=['CANTÓN', 'DENSIDAD_RED_VIAL'],
legend=True,
legend_kwds={
'caption': "Densidad de la red vial",
'orientation': "horizontal"
},
)
# Añadir la red vial al mapa
#red_vial_gdf.explore(
# m=m, # se usa el mapa que se creó en la instrucción anterior
# name='Red vial',
# popup=True
#)
# Agregar un control de capas al mapa
folium.LayerControl().add_to(m)
# Mostrar el mapa interactivo
m
Make this Notebook Trusted to load map: File -> Trust Notebook